Condensation and coexistence in a two-species driven model 
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Condensation transition in two-species driven systems in a ring geometry is studied in the case 
where current- density relation of a domain of particles exhibits two degenerate maxima. It is found 
that the two maximal current phases coexist both in the fluctuating domains of the fluid and in the 
condensate, when it exists. This has a profound effect on the steady state properties of the model. 
In particular, phase separation becomes more favorable, as compared with the case of a single 
maximum in the current-density relation. Moreover, a selection mechanism imposes equal currents 
flowing out of the condensate, resulting in a neutral fluid even when the total number of particles of 
the two species are not equal. In this case the particle imbalance shows up only in the condensate. 
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Many properties of the phase diagram of driven sys- 
tems are knovi^n to be determined by some overall fea- 
tures of the current-density relation. For example, this 
relation serves as a starting point for analyzing models 
of vehicular traffic [1], where it is termed the fundamen- 
tal diagram. It is also a useful tool for analyzing bound- 
ary induced phase transitions in one-dimensional systems 
[2, 3], and stability of shocks [4]. The aim of this pa- 
per is to investigate how these global features affect the 
properties of condensation transitions in driven diffusive 
systems (DDS) on a ring. To this end we analyze in de- 
tail the case where the current-density relation has two 
degenerate maxima. This is found to have far-reaching 
consequences on the emergence of phase separation. It 
results in new features which are not present in the pre- 
viously studied case of a current-density relation with a 
single maximum [5-7]. 

Condensation transitions in one-dimensional DDS have 
been studied in detail in recent years [8]. In particular, it 
was suggested that on a mesoscopic level one can describe 
the dynamics of a broad class of two-species DDS by a 
zero-range process (ZRP) [5, 7]. In this description one 
views the microscopic configuration of the model as a se- 
quence of particle domains, bounded by vacancies. Each 
domain is defined as a stretch of particles of both types. 
Neighboring domains exchange particles through their 
currents. The existence of condensation in these mod- 
els, analogous to Bose-Einstein condensation (BEC), was 
found to be related to the dependence of these currents 
on the length of the domains. A quantitative criterion for 
the existence of a condensation transition in a ring ge- 
ometry was thus suggested. According to this criterion, 
if the asymptotic form of the current for large domains 
of length n behaves as j„ ~ Joo(l + b/n), with 6 > 2, a 
condensation transition takes place at a sufficiently high 
overall particle density. This form of the current implies 
that, at criticality, the steady-state domain-size distribu- 



tion scales as p„ 



for large n. As in the BEC the 



condensed phase is composed of a critical fluid of fluc- 
tuating domains coexisting with a single macroscopically 
large condensate. 

The criterion has previously been applied to models 



where the current-density curve jociv) the bulk of a 
domain exhibits a single maximum [5-7] . Here we apply 
this approach to a model where joo {ri) has two degener- 
ate maxima, and examine the condensation transition in 
cases where, on average, the density within the domains 
lies between the two extremal values corresponding to the 
two maxima. A simple physical picture for the dynamics 
inside a domain is inferred from numerical simulations 
of the model. This picture is substantiated by analyzing 
the properties of a two-species ZRP for modelling the 
collective dynamics of domains. Our main findings are: 
(i) The density in each particle domain (whether a fluid 
or a condensate) is not homogeneous. The two maxi- 
mal current phases coexist within each domain, with a 
sharp interface separating the two. The density profile 
of each of these phases is algebraic, as expected for max- 
imal current phases, (ii) The non-homogeneous density 
profile affects the finite size correction of the current, 
leading to a finite-size correction coefficient B which is 
larger than the expected b of homogeneous domains. For 
example, when the number of particles of both species 
are equal, we find B > 2b. This makes phase separation 
in this model more favorable. Exact solution of the ZRP 
in mean-field geometry and numerical simulations of the 
one-dimensional ZRP support this finding, (iii) In the 
condensed phase the fluid domains are neutral, even in 
systems with non-equal number of particles of the two 
species, leaving the condensate as the only imbalanced 
domain. This is in contrast with the case in which the 
current-density curve has a single maximum, where all 
domains, fluid and condensate, have the same average 
density. 

We now define the model. Consider a one-dimensional 
ring with L sites. Each site i can be either vacant (0) 
or occupied by a positive (-I-) or a negative (— ) particle. 
Positive particles are driven to the right while negative 
particles are driven to the left. In addition to the hard- 
core repulsion, particles are subjected to short-range in- 
teractions through a potential 
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FIG. 1: (a) The current-density relation jca{e,rj) for e = 0.5 
(dashed hne) and e = —0.9 (sohd line), (b) A schematic 
picture of a typical snapshot of the density profile within a 
particle domain. 

where Si = +1 (—1) if site i is occupied by a + (— ) 
particle, and Si = if site i is vacant. The interaction 
parameter e satisfies — 1 < e < 1 to insure positive tran- 
sition rates. The evolution of the model is defined by 
a random-sequential local dynamics, whereby a pair of 
nearest-neighbor sites is selected at random, and parti- 
cles are exchanged with the following rates: 

H — — h with rate 1 -I- AV 

-f ^ + with rate 1 (2) 

— — with rate 1 . 

Here AV is the difference in the potential V between 
the final and the initial configurations. This dynamics 
conserves the number of particles of each species, 7V+ 
and 7V_, or, equivalently, the overall particle densities 
in the system, p± — N±/L. For a given domain, i.e., 
a sequence of positive and negative charges confined by 
vacancies, the relative density r] is simply the fraction 
of positive particles in that domain. This is a fluctuat- 
ing quantity, both in time and from domain to domain. 
Model (2) was studied on a ring geometry for positive e in 
[6, 7]. Here we focus on the negative e region, where the 
current-density relation exhibits two degenerate maxima, 
and study mainly the case p+ = P- • The non-equal den- 
sity case is briefly considered at the end of this Letter, 
and is studied in detail in [9]. 

The region e < was studied in [3] in open systems 
with the purpose of analyzing boundary-induced phase 
transitions. The current-density relation of a domain of 
particles, joo(e, '7), was found to display a single max- 
imum at ?7 — 1/2 for e > —0.8, and two degenerate 
maxima at riH,L = i {l ± [3 - 2((e - l)/e)i/2]i/2} for 
e < —0.8, as depicted in Fig. 1(a). 

We carried out direct numerical simulations of the 
model for e < —0.8. In Fig. 2 we present the domain- 
size distribution and typical configurations for e = —0.9 
at high densities p+ = p_ . This figure suggests the exis- 
tence of a pronounced macroscopic domain. Examining 
the configurations it is evident that the relative density 
within the domains is not homogeneous. Rather, it ex- 
hibits two coexisting regions, corresponding to the two 
maximal-current phases. Indeed, the densities of the two 
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FIG. 2: Domain size distribution for the case e — —0.9. Sim- 
ulation was performed on a system of size L = 5000 with 
7V+ = N- = 1500. Configurations of the model during coars- 
ening towards the steady-state are presented in the inset. 
Positive particles are marked in black, negative particles in 
white, vacancies in gray. A configuration is presented ev- 
ery 500 Monte-Carlo sweeps. The system size is L = 1000, 
N+ = N- = 300, and e = -0.9. 

coexisting phases are equal to riH,L, and the current in 
the system is jooivn) — jooiVL)- This should be com- 
pared with an open system driven with large boundary 
rates, where the system assumes its maximal current and 
a similar coexistence takes place [2]. 

By itself, the appearance of a macroscopic domain in 
numerical simulations of finite systems does not prove 
that condensation takes place, as the presence of such 
a domain could result from a finite size crossover [10]. 
The real question is whether the macroscopic domain 
survives in the thermodynamic limit and becomes a gen- 
uine condensate. To answer this question we use the 
criterion for phase separation, and calculate the finite 
size correction to the current of large domains. When 
the current-density relation has a single maximum, the 
current of a domain of length n takes the asymptotic 
form jn ~ [1 + 6(e, ry)/?!], where ble,ri) is explic- 
itly known [7], and where in all domains 77 is given by 
N^/{N^ + N_). In the present case b must be com- 
puted at the values of the density corresponding to the 
two maxima of the current, rj = rju or 77 = rj^. For 
example we find 6(e = —0.9, 77 = Vh.l) — 1.14. Apply- 
ing the criterion with this value of b would then lead to 
the conclusion that the existence of a macroscopic do- 
main in Fig. 2 is merely a finite-size effect. However, as 
explained below, when a domain is composed of two co- 
existing phases, the real finite size correction coefficient 
which enters the expression of the current is not 6, but 
an enhanced coefficient B > 26, making phase separation 
more favorable. In other words, jn ~ joa ~^ B{e)/n], 
with joo — ioo(e,?7 = Vh,l)- In particular, for e = —0.9 
this yields B > 2, implying that Fig. 2 corresponds to a 
genuine phase separation. 

We first provide numerical evidence that indeed B > 
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FIG. 3: The effective coefficient 6off(e) measured in open sys- 
tems of size L = 400,800, 1600,3200 (from top to bottom). 
The lines correspond to 6(e, rj) (solid line) and 26(e, r;) (dashed 
line), with r) = 1/2 for e > —0.8 and r) = r)H,L for e < —0.8. 

26. It is convenient to calculate tlie finite size correc- 
tion to the current B{e) by simulating an isolated open 
domain of a fixed length n [5]. The coefficient B is 
then extracted by measuring the effective coefficient at 
finite length besin) = n{jn/joo — 1) and extrapolating 
to n — » 00. In Fig. 3 we present bes{n) for various val- 
ues of e and system lengths. It is found that, while for 
e > —0.8 the quantity b^g approaches b{e,r] = 1/2) at 
large n, it is larger than b{e,r] = r)H,L) by a factor > 2 
for e < —0.8. We note that higher order corrections be- 
come significant as one approaches e = —0.8, where the 
leading finite-size correction, 6(e,r7 — l/2)/n, vanishes. 

We now present a physical explanation of these ob- 
servations. For a domain of length n each of the two 
coexisting phases occupies on average only a length n/2. 
This effectively reduces the length of the domain by a fac- 
tor 2, and thus the finite size correction is expected to be 
about ~ 26/n rather than b/n. Quantifying this intuitive 
picture leads to an estimate of the enhancement factor 
A = B{e)/b{e,r]H,L)- We analyze the current emitted 
from a domain of length n in the fluid, composed of two 
coexisting maximal-current phases. A schematic density 
profile in such a domain is given in Fig. 1(b). At the 
left side of the domain a fraction x of its length is oc- 
cupied by a phase with high bulk density rjH, while the 
remaining right side is occupied by the other maximal- 
current phase, with low density rjL- The position of the 
interface fluctuates in time around the midpoint, i.e., on 
average, {x) = 1/2. Numerical simulations strongly sug- 
gest that the position x varies on time scales which are 
much longer than the equilibration time of the local den- 
sity within each phase [9] . We thus consider the dynamics 
of the system on time scales which are short enough such 
that the position of the interface x and the size of the 
domain n are fixed. On these time scales the currents of, 
say, positive particles jnix) and jL{x) in the high and 
low density phases, respectively, are given by 

Jh{x) = Joo (l + , Mx) = Joo (l + ;^(^) . 

(3) 



Thus, as a result of the flow of particles through the 
domain, the interface moves with a velocity v, such that 
3h{x) — jhix) = V {r}H — Vl) ■ The outflow of particles 
from the domain is therefore given by 

jnix) - vrjH = jhix) - VT]L = joo [1 + A{x)b/n] , (4) 

where, using the expressions above, one has 

A(-) = ^— (t^-^) • (5) 

Vl-Vr \^-x X J 

On longer time scales where the position of the interface 
X fluctuates one has to average (4) in order to get the 
current emitted from the domain, leading to A = {A{x)). 
If the fluctuations in the position of the interface do not 
scale with the domain size, then (l/x) = l/{x) = 2, and 
hence A = 2. On the other hand, if these fluctuations 
scale like the domain length, then (l/x) > l/{x), and 
A> 2. In the following we explore this question in more 
detail. Our analysis suggests that indeed the width of 
the interface scales with the domain length leading to 
B > 26. 

Motivated by the discussion presented above, we in- 
troduce a two-species ZRP which captures the main fea- 
tures of the collective dynamics of the evolving domains. 
Consider a ring of M boxes, where box i contains rij 
particles, ki of which are positive and k are negative: 
rii = ki + li. The dynamics of the model is such that 
a box i is chosen at random and a positive charge is 
moved to its right neighboring box with rate Uk,i and a 
negative charge moves to its left neighboring box with 
rate Vk.i- In this model a box represents a generic do- 
main of the original DDS, and the rates u and v corre- 
spond to the outflow of particles from this domain, as 
found in (4). We thus take Ukj = 1 + A{x)b/n and 
Vk,i = 1 + ^(1 — x)b/n. The variable x relates to the 
relative density r] by xr]L -|- (1 — x)r]ii = r] = k/n. In 
what follows we analyze for simplicity the case Jji = 1 
and r]R = 0, which yields 

6 6 
Uk,i = l + j, Vk,i = 1 + - . (6) 

With this choice of rates the steady state of the model 
is not a product measure [11], implying that no explicit 
description of the stationary state is known. We first 
consider the model in the mean-field geometry, where all 
sites are connected. Wc denote by fk.i the probability 
for a site to be occupied by k positive particles and I 
negative particles. In the thermodynamic limit fk,i obeys 
the evolution equation 

r- = Uk+lJ ,tk+l,l + Vk.l + 1 Jk,l+1 

at 

+ u fk-id{l - Sk,o) + V ,fk,i-i{l - Si^n) (7) 
- - h,o) + Vk,i{l - Sift) +u + v] fk,i , 
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FIG. 4: The decay exponent a for different values of h, as 
obtained from numerical integration of (9). The line is given 
by the large 6 asymptotic form 26 + %/2. 
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FIG. 5: The relative density 77 of positive particles within 
a domain, for a system of size 10"* with e = -0.9. Here 
N+ = 5000 and iV_ = 3000. 



where u = J^kj'^kdfkd and v = Y.k,i'"k,ifkj are the 
(+) and (— ) currents, respectively. In the continuum 
limit the steady state distribution satisfies the following 
equation at criticality {u = v= 1) 



^ dp 



1 dfk,i _^ 1 dfkd 



I dk k dl 



. 



(8) 



Moving to polar coordinates, and assuming the scaling 
solution f{r,6) = r~'°'g{9), we find an equation for the 
angular function g{6) 



d^2 



+ U- 



26 
sin 261 



ag{e)=Q, 



(9) 



with the boundary conditions 5(0) = 5(7r/2) = 0. The 

determination of the decay exponent a is obtained by im- 
posing the boundary conditions. This is the quantization 
condition for this Schrodinger equation. Except for spe- 
cial values of b where a can be determined exactly (e.g. 
a = 3 for 6 = 2/3), the value of a as a function of b is ob- 
tained by integrating (9) numerically (Fig. 4). The large 
b asymptotic form obtained by the WKB approximation, 
26-F V2, agrees very well with these results down to 
small b. From the predicted form of the solution /(r, 6) 
we deduce that the domain size distribution p,i, with 
n = A; -t- /, scales as pn ~ n~^°'~^\ On the other hand, 
the rate out of a domain of size n, jn = {iJ'k,i)k+l=n, is of 
the form jn ^ 1 + B/n, and wc conclude that B = a — 1. 
Numerical integration of the temporal eqs. (7) gives a de- 
cay exponent a in perfect agreement with the predicted 
value of the continuum limit. This analysis shows that 
B > 2b, supporting the physical picture presented above. 
This calculation was carried out within the mean-field 
geometry and should not yield the exact values of B of 
the one-dimensional model. However, numerical simula- 
tions of the latter indicate that B is well approximated 
by the mean-field result [9]. Coming back to the DDS, 
the results above suggest that the position of the inter- 
face inside a domain should scale with the domain length. 
This has been verified by numerical simulations [9]. 

So far we analyzed neutral system, where p+ = P-- 
We now consider the case of non-equal densities. While 
this case will be studied in detail elsewhere [9], here 
we only mention a striking result: as long as r]L < 



N+/{N+ + N-) < r]H, all domains which reside in the 
fluid are equally populated with positive and negative 
particles. The excess number of particles of the major- 
ity species reside in the condensate. This behavior is a 
result of the fact that the two currents emitted from a 
domain of length n are equal to the maximal current, up 
to corrections of order 1/n. The condensate is therefore 
stationary in the thermodynamic limit even when the 
densities are not equal. The condensate thus emits equal 
currents of (-I-) and (— ) particles. Hence domains in the 
fluid cannot experience the fact that the overall densities 
in the system are not equal. Fig. 5 depicts the average 
relative density 77 of positive particles in domains of var- 
ious sizes, as measm'cd in a large system for e = —0.9. It 
is readily seen that on average domains in the fluid are 
neutral, whereas the relative density in the macroscopic 
domain compensates for the excess of positive particles. 
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